function [ ] = plotFitHistEmissionsAll( plant, ctrM, scenarios_toplot, ...
    filestub, titlestub, inlogs, optfig )
% Plot total costs under the emissions market and under the command and
%   control scenario
%
%   INPUTS:
%     plant - Table of emissions and other outcomes
%     ctrM  - Counterfactual object with data simulated from estimated
%             parameters
%     scenarios_toplot - Table of scenarios to plot
%     filestub - Root of filename in which to print figure
%     titlestub - Root of title to give the figure
%     inlogs - Whether to plot emissions histogram in logs
%     optfig - Figure options

% Pick out only the scenarios passed as arguments to plot
scenarios_toplot = renamevars(scenarios_toplot,'cap','id_level');
ctrTable = ctrM.output_itk;
ctrTable = innerjoin(ctrTable,scenarios_toplot,...
      'Keys',{'id_period','id_level'});

tabulate(ctrTable.id_period);

% Get axis limits across all periods
Eitall = [ plant.Eit; ctrTable.Eit ];
Elim = prctile(log(Eitall),[0.1 99.9]);
Elim(1) = floor(Elim(1));
Elim(2) = ceil(Elim(2));

% Plot histogram for each period
for i = 1:size(scenarios_toplot,1)
    
    % Get parameters describing the scenario: period, cap
    t   = scenarios_toplot.id_period(i);
    cap = scenarios_toplot.id_level(i);
    
    % In data
    keep = (plant.id_period == t);
    Eit  = plant.Eit(keep);
    apcd_max=false;
    filename = [ filestub sprintf('data_%02.0f',t) ];
    title    = [ titlestub sprintf(': data, period %2.0f',t) ];
    title = "";
    plotFitHistEmissions( Eit, apcd_max, filename, t, title, inlogs, optfig, Elim, false);
      
    % In simulation (counterfactual)
    keep = (ctrTable.id_period == t & ctrTable.id_level == cap);
    Eit  = ctrTable.Eit( keep );
    filename = [ filestub sprintf('simul_%02.0f',t) ];
    title    = [ titlestub sprintf(': counterfactual, period %2.0f',t) ];
    title = "";
    apcd_flag = false;
    % To have the apcd type breakdown for step functions uncomment the below
    %apcd_flag = ctrM.scenario.costs=="step"; 
    if apcd_flag
        spare = ctrTable;
        spare.apcd_max(spare.minPeriodEmissions==spare.E_abate) = -1+0*spare.apcd_max(spare.minPeriodEmissions==spare.E_abate);
        apcd_max = spare.apcd_max(keep);
    end
    plotFitHistEmissions( Eit, apcd_max, filename, t, title, inlogs, optfig, Elim, apcd_flag);
    
end

end